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I.  MOTT'S  DISTRIBUTION 


During  World  War  II  in  England,  Mott  and  Linfoot1  tried  represent¬ 
ing  the  number  of  fragments  produced  by  the  detonation  of  a  shell 
or  bomb  by  the  distribution 

dN  =  [N^y)]  (m/y)"2/3  exp  [-  (m/y)1/3j  dm 

=  Nt  exp  [-  (m/y)1//3]  d  [(m/y)1>/3]  (1) 

where  dN  is  the  number  of  fragments  in  the  mass  interval  from  m  to 
m+dm.  Here  NT  is  total  number  of  fragments  as  can  easily  be  seen 
by  letting  t=(m/y)  J'/  ^  and  integrating  Equation  (1)  from  0  to 
namely. 


oo  /*  CO 

dN  =  Nt/  e_t  dt  =  Nt.  (2) 

T*/  o  T 

The  parameter  y  can  be  related  to  the  expected  value  of  the 
distribution,  m,  as  follows: 


m 


-r 


-j 

•'a 


m  (.dN/N  J  =/  yt  (e-tdt)  =  y  (31).  (3) 

o  1  Jo 

They  used  an  experimental  statistic,  namely,  the  observed  average 
mass  of  the  fragments  collected,  as  an  estimator  of  the  expected 
value  in  to  obtain  an  estimate  of  the  population  parameter,  y  =  m.,,/3!), 
where  m„r  =  M^/N,,,  with  M^,  being  the  total  fragment  mass  collected! 

This  gave  rough  agreement  with  experimental  data  obtained  for  a  3.7  inch 
shell  for  example,  at  least  for  the  smaller  fragments.  For  the  larger 
fragments  (which  often  showed  part  of  the  original  shell  inner  and 
outer  walls  as  surfaces)  they  observed  that  a  better  description 
could  be  obtained  by  using  1/2  instead  of  1/3  in  Equation  (1).  This 
suggested  to  them  a  dominance  of  three-dimensional  fracture  in  the 
production  of  small  fragments  and  a  dominance  of  two-dimensional 
fracture  for  larger  fragments,  leading  to  a  different  estimate  for  y, 
namely,  y  =  mA„/  (2!).  In  applications  Mott  never  actually  split  a 
population  into  two  groups,  using  exponents  1/3  and  1/2,  but 
used  one  or  the  other.  He  remarked  that  similar  distribution  laws 
had  already  been  applied  to  the  crushing  of  rock,  and  cited  an 


1  N.F.  Mott  and  E.  H.  Linfoot,  "A  Theory  of  Fragmentation,"  A.  C. 
3348,  Jan  1943. 
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American  reference2. 

In  a  later  report,  Mott3  suggested  that  for  the  exponent  1/2  at 
least,  the  parameter  y  could  be  estimated  from  the  shell  wall 
thickness,  T,  and  its  inner  diameter,  D,  by  the  formula 

M1/2  =  C  T5/6  D1/3  (1  +  T/D)  (4) 

where  the  empirical  constant  C  is  larger  for  explosives  which  impart 
a  lower  launclji  .^elocity  to  the  fragments.  He  did  not  give  a  similar 
formula  for  y  ,  the  form  found  for  this  parameter  in  Equation  (1) . 
Since  real  shells  have  variable  T  and  D  and  C  along  their  length,  it 
is  clear  that  Mott  was  here  considering  the  idealized  case  of  a 
right  circular  cylinder.  Later4  be  considered  this  idealized  case 
of  a  cylinder  formed  from  stacked  metal  rings  in  more  detail.  Here 
fracture  perpendicular  to  the  shell  axis  is  pre-determined,  and 
only  ring  fracture  need  be  considered.  This  is  similar  to  the 
natural  fragmentation  of  a  real  warhead  if  we  look  on  the  shell  as 
made  up  of  stacked  rings  of  different  size  and  shape. 

Shortly  after,  a  co-worker  of  Mott  named  Ursell5  suggested  that  the 
one-dimensional  fracture  of  a  rod  ought  to  be  given  by  a  Poisson  dis¬ 
tribution,  and  might  be  related  to  warhead  fragmentation.  Still  later 
in  the  United,  States,  Thomas6  pointed  out  formulas  of  the  type 
exp  [- (m/y) 1/( *]  where  £  =  1,  2  or  3  for  one-,  two-  or  three-dimensional 

fracture  are  merely  probability  distributions  more  or  less  suited  to 
describing  particular  fragment  populations.  It  is  not  very  helpful 
to  require  a  physical  model  which  envisions  simultaneous  formation  of 
fragments  by  means  of  planes  traversing  volumes,  lines  traversing 
planes  or  points  dividing  lines.  More  realistically,  all  of  these 


2  Lienau,  J.  Franklin  Inst.  p485  ("1935). 

3  N.  F.  Mott,  Fragmentation  of  H.  E.  Shells;  a  Theoretical  Formula 
for  the  Distribution  of  Weights  of  Fragments,  A.C.  3642,  1943. 

4  N.  F.  Mott,  "Fragmentation  of  Shell  Cases,"  Proc.  Roy.  Soc. 
(Lond.)  189,  300  (1947). 


5  H.  D.  Ursell,  "Fragmentation  Data  and  Theories  of  Fragmentation," 
A.  C.  3817,  1943. 

6  L.  C.  Thomas,  "Comments  on  Mott's  Theory  of  the  Fragmentation  of 
Shells  and  Bombs,"  BRL  R398,  Sept  1943,  (AD  #36152) 
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processes  might  be  going  on  simultaneously  and  sequentially,  with 
smaller  fragments  being  formed  later  from  larger  fragments  already 
formed.  In  spite  of  this,  Thomas  used  l  =  2  in  his  own  applications 
to  U.  S.  munitions,  possibly  because  of  the  calculational  difficulties 
involved  with  the  use  of  non- integer  Jt  at  the  time.  Thomas. seems  to 
have  been  the  first  to  describe  the  fragmentation  of  real  warheads 
by  applying  Mott's  formula  to  individual  rings  of  different  wall 
thickness  and  diameter  formed  by  mentally  slicing  a  shell  perpendicular 
to  its  axis7.  This  practice  is  still  in  use  today8. 

G.  I.  Taylor  also  considered  the  explosive  fragmentation  of  metal 
rings  with  radial  cracks  starting  on  the  outside  of  the  shell  case 
and  propagating  inward  under  the  combined  influence  of  tensile  and 
compressive  forces9.  He  also  pointed  out  that  cracks  should  propagate 
at  about  45°  to  the  circumferential  and  radial  directions,  further 
complicating  fragment  size  and  shape  distributions.  More  complete 
treatments  such  as  that  of  Nadai10  indicate  that  a  system  of  logarithmic 
spiral  cracks  should  develop  in  the  wall  in  the  simple  case  of 
infinitely  long  cylinders  of  uniform  wall  thickness  uniformly  stressed 
from  the  interior.  Of  course,  for  real  warheads  the  crack  systems  will 
be  much  more  complicated.  A  further  complicating  factor  is  the 
existence  of  shock  waves  which  reverberate  in  the  shell  wall  as  it 
expands 1 1 

In  summary,  the  physics  of  real  warhead  fragmentation  is  so 
complicated  that  simplified  models  are  not  likely  to  be  much  help  in 
predicting  fragment  size  distributions.  When  confronted  with  problems 
of  extreme  complexity,  physicists  generally  invoke  some  form  of 
probability  description  as  in  statistical  mechanics.  In  what  follows 
we  will  try  to  point  out  that  Mott's  general  procedure  can  be  given 
a  rational  basis  in  statistical  theory  and  can  be  improved  somewhat 


7  L.  H.  Thomas,  "Analysis  of  the  Distribution  in  Mass,  in  Speed, 
and  Direction  of  Motion,  of  the  Fragments  of  the  M71  (90mm)  A.  A. 
Shell,  when  Filled  with  TNT  and  when  Filled  with  Ednatol,"  BRL 
R434,  Dec  1943. 

8  Glenn  Rand ers-Pehr son,  R.  R.  Karpp,  C.  E.  Anderson,  Jr.  and  H.  J . 
Blische,  "Shortfrag  Users  Guide,"  ARBRL-MR-03007,  Mar  1980. 

9  G.  I.  Taylor,  "The  Fragmentation  of  Tubular  Bombs,"  in  TTie 
Scientific  Papers  of  Sir  Geoffrey  Ingram  Taylor,  Ed.  by  G.  K. 

Batchelor  (Cambridge:  The  University  Press,  1963)  v.3,  p387. 

10A.  Nadai,  Theory  of  Flow  and  Fracture  of  Solids,  (N.Y.:  McGraw- 
Hill,  1950)  p539 . 

11F.  E.  Allison  and  J.  T.  Schriempf,  "Explosively  Loaded  Metallic 
Cylinders,  II,"  J.  Appl.  Phys.  31,  846  (1960). 
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once  it  is  realized  that  we  are  using  a  general  theory  of  random 
breakup  rather  than  particular  physical  models.  We  will  also  point 
out  the  generality  of  the  description  by  showing  that  it  can  be 
applied  to  controlled  as  well  as  to  natural  fragmentation.  We  will 
not  attempt  to  link  microscopic  fracture  theories  to  continuum  fracture 
mechanics  by  constructing  distribution  functions  which  represent 
particular  types  of  flaw  activation  rates  as  others  have  done12. 
Instead,  we  will  postulate  an  overall  defect  activation  rate  and  pursue 
the  statistical  consequences  of  such  an  assumption.  In  particular, 
we  will  inquire  how  general  and  simple  we  can  keep  the  form  of  our 
distribution  function  without  losing  its  ability  to  represent 
experimental  data  well  enough  for  practical  purposes. 


II.  A  GENERAL  DISTRIBUTION 
A.  Derivation  of  the  Distribution. 

Consider  a  solid  body  of  volume  y  and  any  shape-  Real  solids 
generally  contain  many  kinds  of  defects  which  can  act  as  weak  points 
when  a  stress  is  applied  and  serve  to  initiate  cracks  and  fractures. 

Let  us  mentally  divide  this  body  into  k  elementary  volumes  of 
size  e  =  y/k,  choosing  k  large  enough  so  that  on  average  each 
elementary  volume  e  contains  one  defect  or  incipient  break.  Let 
r  be  the  average  volume  rate  of  defect  activation  under  stress  so 
that  re  =  ry/k  is  the  probability  of  finding  at  least  one  such 
activated  defect  in  e.  Then  the  probability  of  finding  no  such 
defects  in  e  is  1-re.  Since  the  distribution  of  defects  and  their 
activation  under  stress  can  be  considered  random  events,  the 
probability  of  observing  exactly  s  such  activations  in  k  trials  is 
given  by  the  binomial  distribution 

*  s:m:  (F)S  (l-^)k’S  <5> 

for  s  =  0,l,2...k  with  B  =  0  otherwise.  If  s=k  for  example,  then 
all  defects  would  be  changed  into  breaks  and  the  body  would  be 
subdivided  as  much  as  it  could  be  by  the  defect  mechanism.  Since 
the  number  of  defects  in  real  solids  of  interest  is  very  large,  it 
should  be  adequate  to  consider  the  limit  as  k  -*•  “  in  such  a  way  that  the 
number  of  breaks  or  activated  defects,  ry,  is  large  but  far  from 
infinite. 


12D.  R.  Curran,  L.  Seaman  and  D.  A.  Shockey,  "Dynamic  Failure  in 
Solids,"  Physics  Today,  Jan  1977,  p.  46. 
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In  particular 


ry/k  <<  ry  <<  k 


(6) 


is  the  condition  of  interest.  In  this  case  from  Equation  (5) 


p(s/ry) 


(7) 


where  p(s/ry)  is  the  Poisson  distribution  with  parameter  ry.  The 
first  limit  in  Equation  (7)  is  seen  to  be  unity  if  we  divide  numerator 
and  denominator  by  (k-s) ! ,  so 


k(k-l)  (k-2).. .  (k-s+1) 
lim  ,  s 

k~>  k 


lim 

k-*» 


*)(>-  r}-(l-tTy)- 


1.  (8) 


The  third  limit  in  Equation  (7)  is  obviously  unity.  If  we  let  u  =  ry 
and  z  =  -u/k,  then  we  see  that  the  second  limit  in  Equation  (7)  is  just 
the  (-u)  power  of  the  limit  which  defines  the  base  of  the  natural 
logarithm,  namely 


(9) 


Of  course  the  Poisson  distribution  meets  the  requirement  that 


])(:-  /u) 

s=o 


s=o 


oo 


s=o 


Let  F(u)  be  the  probability  that  at  least  one  break  (activated 
defect)  will  occur  under  a  given  stress.  Then  l-F(u)  is  the  probability 
that  no  break  will  occur.  From  Equation  (7),  the  probability  of  no 
breaks  (s=o)  is  l-F(u)  =  e-u,  or 


F(u)  =  1-e  U  P (s/u)  -  e“U  jr  e'U  (11) 

s=o  s=l 
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which  exhibits  the  meaning  of  F(u)  as  the  probability  of  observing  at 
least  one  (one  or  more)  break.  We  note  in  Equation  (11)  that  F(u) 
has  been  written  both  as  a  continuous  function  of  the  variable  u  and 
as  a  sum  over  a  discrete  frequency  distribution  in  which  u  plays 
the  role  of  a  parameter.  This  is  a  particular  case  of  a  more  general 
relation  as  we  shall  see  below.  Of  course  F(u)  can  also  be  written 
as  an  integral  over  a  continuous  frequency  distribution  as  follows: 


where 


=  J  f  (t)dt 

rv 


(12) 


dF  =  e_t  dt  =  f (t)dt. 


(13) 


Here  f(t)  is  a  frequency  distribution  and  F(u)  is  a  cumulative 
distribution,  while  1-  F(u)  is  a  complementary  cumulative  distribution. 
Equations  (11)  and  (12)  are  both  particular  cases  of  more  general 
relations,  namely. 


F  (u,  c) 


-£  £  ^  Xfef  ■  fU  rfc 

s=c  ^  o 


c-1 


dt.  (14) 


Equations  (11)  and  (12)  are  obtained  for  the  case  c=l  in  Equation  (14) 
since  the  complete  gamma  function  of  unity  is  unity,  r ( 1)  =  Oi  =  1. 

The  complementary  form  of  Equation  (14)  is 


c-1 

1-F(u,c) 

s=o 


u 
s ! 


e  U  =  Q(2u/2c) 


r(c,u) 

r(c) 


/*-  l 

Vu  r(c) 


c-1  - 1 

tC  e  ldt  (15) 


and  is  discussed,  for  example,  in  the  Bureau  of  Standards  Handbook 
of  Mathematical  Functions13  where  tables  of  these  functions  are 
also  given.  Equation  (14)  relates  partial  sums  over  the  Poisson 
distribution  to  the  chi-square  distribution,  p(2u/2c)  with 

_  SL=1  ' 

13  Handbook  of  Mathematical  Functions,  NBS  Applied  Mathematics  Series, 
number  55,  Nov.  1964,  Ed.  by  M.  Abramowitz  and  I.  A.  Stegun,  section  26, 
especially  26.4.19,  26.4.21  and  26.4.2. 
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where  2c  is  called  the  number  of  degrees  of  freedom.  As  we  see, 
these  sums  are  also  related  to  the  incomplete  and  complete  gamma 
functions,  y(c,u)  and  r(c).  The  integrand  in  Equation  (14)  is  the 
gamma  frequency  distribution  with  scale  factor  unity  and  is  a 
special  case  of  the  Pearson  type  III  distribution: 


1  c-1  -t 

hi  er(c)  e 


where 

*  =  O  -  x0)/e] 


(17) 


(18) 


for  y  <  y  <  oo  and  scale  factor  3  [NBS  Handbook,  p930]  .  An  even 
more  general  form  of  Equation  (18)  is 


t  =  L(y  -  y0)/e] 


1A 


(19) 


where  0<2.<°°.  This  gives  us  Weibull's  frequency  distribution 
dF  =  W (t)  dt  =  e_t  dt  =  exp  ~  d  I j 


dy  =  W(y)  dy  (20) 


which  contains  the  exponent  1/i.  as  well  as  the  scale  factor  8  and 
the  cutoff  value  y  for  a  random  variable  y  over  the  range  y  <y<oo. 
Weibull11*  has  appl?ed  this  distribution  to  a  great  variety  o? 
phenomena,  social  as  well  as  physical,  chemical  and  biological,  as 
have  others  after  him.  Gnedenko15  had  shown  previously  that  Equation  (20) 


14  W,  Weibull,  "A  Statistical  Distribution  Function  of  Wide 
Applicability,"  J.  Appl.  Mech.,  Sep  1951,  p293. 

15  B.  V.  Gnedenko,  "Limit  Theorems  for  the  Maximum  Term  of  a 
Variational  Series,"  Doklady  Akad.  Nank,  USSR32,  1941. 
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is  the  third  asymptotic  distribution  of  smallest  values16. 

We  note  that  Mott's  distribution  in  Equation  (1)  is  a  particular  case 
of  Equation  (20)  for  Jl  =  3,  y  =  m,  yQ  =  m  =  o,  3  =  m  and  d(N/T  )  = 
W(t)dt  =  dF.  Mott's  formula  with  i  =  2  i?  another  special  case1 
of  an  integer  l  value.  Since  Equation  (20)  is  a  particular  case 
of  Equation  (14)  with  c  =  1  and  t  given  by  Equation  (19),  so  is 
Mott's  formula.  A  more  general  form  of  Weibull's  distribution 
which  may  be  used  in  Equation  (14)  with  any  allowable  c-value  is 
the  general  gamma  distribution: 


E(t)  dt 


g(y)  dy 


/  \l/£ 

/  „  \i/* 

X 

1 

X 

o 

A 

I"  o) 

\  3  / 

a 

V  e  / 

(V-F 

dy 

(21) 


Here  y  is  a  general  random  variable,  but  in  fragmentation  applications 
we  take  it  to  be  a  volume,  y  =  m/p,  where  p  is  the  density  of  the 
shell  case  and  m  is  mass.  With  the  scale  factor  3  =  p/p,  c  =  1  and 
l  =  2  or  3,  Equation  (21)  gives  Mott's  formulas.  Since  c  =  1  in 
Mott's  formulas,  we  see  that  he  is  calculating  the  probability  of  at 
least  one  break  occurring.  We  have  also  seen  that  it  is  a  special 
form  of  general  probability  distributions  and  can  be  applied  to  many 
things  besides  fragmentation.  As  we  mentioned  above,  Thomas  pointed 
out  that  Mott's  formulas  are  not  necessarily  connected  to  any 
simple  model  of  the  fragmentation  process.  More  generally,  we  now 
see  that  they  are  not  necessarily  limited  to  fragmentation  at  all. 

Let  us  return  to  Equation  (14)  and  display  some  particular 
examples  by  way  of  illustration.  If^the  parameter  u  =  1,  the 
Poisson  frequency  function  is  p  =  e”  /(s.'),  which  is  plotted  in 
Figure  1.  In  this  case  the  chance  of  obtaining  one  break  (s  =  1)  is 


16  E.  J.  Gumbel,  "Statistical  Theory  of  Extreme  Values  (Main  Results)," 
c.6  in.  in  Contributions  to  Order  Statistics,  Ed.  by  A.  E.  Sarken 
and  B.  G.  Greenberg  (N.Y.):  John  Wiley  and  Sons,  Inc,  1962. 
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the  same  as  the  chance  of  obtaining  no  breaks  (s=0) ,  an  instance  of  the 
double  mode  or  maximum  of  the  Poisson  formula.  As  is  well  known,  the 
expected  value  of  s  is  equal  to  the  parameter  u  (as  is^  the  variance) 
for  the  Poisson  distribution.  That  is,  in  Figure  1,  s  =  u  =^1. 

Since  the  chance  of  obtaining  no  breaks  is  large,  namely,  e“  =  0.368, 
the  chance  of  obtaining  at  least  one  break  (  one  or  more  breaks)  is 
only  0.632.  The  chances  of  obtaining  s  =  2,3,4....  breaks  become 
smaller  as  s  increases  and  the  chance  of  obtaining  5  or  more  breaks 
is  quite  small.  From  Equation  (14)  it  is 

OO  4 

F (1 , 5)  ^2  e_1/Cs!)  =  1^2  e_1/(s!)  =  1"e 
s=5  s=0 

=  y  (5,l)/r(5)  -f  j,  t4  e-t  dt  =  — I?---  =  .00366  .  (22) 


I;f  the  parameter  u  is  larger,  say  u  =  5,  the  Poisson  frequency 
is  e"“  5s/ (sj)  which  is  plotted  in  Figure  2.  Now  thg  chance  of 
obtaining  no  breaks  at  all  is  quite  small,  namely  e  =  0.0067,  but 
is  not  zero.  The  expected  number  of  breaks  is  s  =  u  =  5  while  the 
probability  of  obtaining  at  least  one  break  is  now  much  larger, 
namely,  from  Equation  (14) 

00  s  r 5 

F (5, i)  =y^  |t-  e‘5  =  Y  (1,5)  =  I  e_t  dt  =  1-e"5  =  0.9933  (23) 

”  J  o 

while  the  probability  of  obtaining  at  least  the  expected  number 
of  breaks  is 

OO 

f(5, 5)  =^2  It  e"5  =  HW  =  h:  I " &  t  dt  =  °-5595  •  (24) 

s=5  J  o 


In  short,  an  increase  in  the  parameter  u  reduces  the  chance  of  no 
breaks  at  all  and  shifts  the  distribution  to  the  right,  since  the 
expected  number  of  breaks  is  always  s  =  u.  A  decrease  in  u  has  the 
opposite  effect.  If  u<l,  s  is  fractional,  and  the  probability  of 
obtaining  no  breaks  at  all  is  the  mode  or  most  likely  event. 
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Figure  2.  Poisson  Frequency  Distribution  with  Parameter  p  =  5. 


In  the  limit  as  u  +  o,  the  Poisson  distribution  for  s  =  o  is  a 
spike  function  equal  to  unity  at  s  =  o  and  zero  for  all  other  s. 

Of  course,  F(0,1)  =  1-e  =  0,  that  is,  the  chance  of  at  least  one 
break  vanishes.  In  other  words  for  a  perfect  crystal  for  which 
u  =  ry  =  0  because  the  rate  of  occurrence  of  defects  r  =  0,  we  are 
sure  that  no  break  will  occur  by  a^defect,  mechanism.  In  the  limit 
as  u  -*■  00  ,  the  poisson  function  (u^ / s.') / eu  is  indeterminate. 
However,  using  L'  Hospital's  rule,  s  differentiations  of  the 
numerator  reduces  it  to  (s!)/(s!)  or  unity,  while  the  denominator 
remains  the  same,  so  that  each  term  vanishes.  Actually,  u  is  not 
allowed  to  increase  without  limit  as  was  pointed  out  in  Equation  (6) 
where  u  =  ry  <  <  k.  Even  though  k  is  allowed  to  increase  without 
limit,  u  must  remain  finite.  'Phis  requirement  agrees  with  our 
interpretation  of  the  dimensionless  number  u  as  the  finite  number 
of  defects  expected  to  be  activated  under  a  given  stress.  It  also 
agrees  with  an  interpretation  of  u  as  the  ratio  of  the  stress 
energy  applied  per  unit  volume  to  the  work  per  unit  volume  required 
to  fracture  the  body,  namely. 


applied  energy/ volume 

at; 


(25) 


where  a  and  e  are  the  stress  and  strain  at  fracture  characteristic 
of  the  material.  In  the  case  of  projectiles  striking  a  target  at 
ordnance  speeds,  u  will  not  be  much  greater  than  unity  and  the 
projectile  will  fracture  into  several  pieces17.  In  this  case 
discrete  Poisson  statistics  are  appropriate.  For  hypervelocity 
impacts,  a  projectile  will  shatter  into  many  fragments,  as  also 
happens  when  explosive  warhead  cases  are  shattered  at  detonation. 

In  such  cases  a  continuous  frequency  distribution  is  more  convenient 
for  describing  the  resulting  fragment  population.  Equation  (14) 
shows  how  these  distributions  are  connected.  When  u  =  1  in  Equation 
(25),  the  applied  force  is  matched  by  the  strength  of  the  body  and 
on  average  we  expect  only  one  break  to  occur,  although  it  is  equally 
likely  that  the  body  remain  intact.  When  u  is  much  less  than  unity 
we  are  in  a  regime  of  slow  crack  development  and  failure  which  may 
take  months  or  years  of  stressing.  For  explosive  ordnance  we  expect 
u  much  greater  than  one  with  fragmentation  times  in  the  millisecond 
range.  But  in  all  cases  u  will  be  finite  and  the  upper  limit  of 
the  integral  in  Equation  (14)  will  not  be  infinite,  although  in  some 
cases  of  interest  it  may  be  effectively  infinite  to  a  sufficient 
approximation. 


17  J.  Dehn,  "The  Particle  Dynamics  of  Target  Penetration," 
ARBRL-TR-02188,  Sept  1979. 
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If  we  use  y  =  m/p  and  3  =  p/p  in  Equation  (19)  we  obtain 


and 


(26) 


(27) 


so 


m  <m<m  =  m  +  p  u  (28) 

o  u  o 


where  m  is  the  upper  mass  limit  of  a  finite  body.  For  u  greater 
than  teK,  we  see  that  the  probability  of  at  least  one  break 
occurring,  F(u)  =  l-eu  is  very  close  to  unity.  Since  the  number 
of  fragments  expected  is  one  larger  than  the  number  of  breaks,  they 
are  approximately  equal  for  a  large  number  of  breaks. 


For  u  =  1  in  Equation  (27),  p  =  (m^  -  mQ)  for  any  A.  In  general, 

for  u  —  1,  p  —  (m  -  m)  for  any  A  (0<A<°°)  .  Here  p  is  relatable 
<  >  u  o 

to  an  expected  average  mass  and  smaller  p  is  associated  with  larger  u, 
that  is  with  greater  applied  stress,  a  weaker  body  and  more  breaks 
or  fragments  expected.  The  terms  for  small  s  in  Equation  (14) will 
be  very  small  for  large  u,  that  is,  the  occurrence  of  only  a  few 
breaks  is  very  unlikely.  Most  of  the  contribution  to  either  the 
summation  or  integration  in  Equation  (14)  will  come  from  the 
midrange  near  s  =  u  or  m  near  p.  For  example  in  Equation  (14)  we  can 
take  c  =  0  or  1  and  obtain  the  expected  number  of  breaks 


■£  4- «)£•(*■£)■  "'■££■ 


-u  u 
=  ue  e  =  u 


(29) 


s=o 


S=1 


3=o 


since  the  tejm  for  s  =  0  makes  no  contribution.  Similarly,  using 
m  =  ra0  +  u  t  from  Equation  (26)  we  can  find  the  expected  values  of 
the  mass  (letting  c  =  1) 
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(30) 


mQ  Y  (l,u)  +  y  y  (1+A,u) 


If  u  is  effectively  infinite,  then  the  incomplete  gamma  function,  y, 
is  approximately  equal  to  the  complete  gamma  function,  r,  and  Equation 
(30)  becomes 


m  =  mn  r(l)  +  y  r  (1+A)  =  m^  +  y  r  (1+A)  .  (31) 


If,  in  addition,  the  lower  mass  limit,  m  ,  is  effectively  zero  and  A 
is  an  integer,  Equation  (31)  becomes  ° 


m  =  y  (A')  (32) 


since  r  (1+A)  =  AJ  for  integer  A.  Mott's  formulas  consider  only 
the  values  A  =  2  or  3.  In  general,  A  need  not  be  integer,  mQ  need 
not  be  zero  and  m^  (and  so  u)  is  less  than  infinite. 

In  Figure  3  we  plot  the  incomplete  gamma  function  y(1+A,u)  versus 
u  for  A=0,l,2,3  and  4.  For  non- integer  A,  the  curves  lie  between 
those  shown.  It  is  clear  that  for  smaller  A  values  the  approximation 
y(l+A,u)  **  r(l+A)  is  quite  good  for  smaller  u  than  for  larger  A  values. 
For  0<A<3,  it  is  very  good  for  u-ilO.  This  was  mentioned  above  in 
another  way  when  we  observed  that  F(u)=l-e”  1  for  u>10. 

B.  Description  of  Fragment  Populations 

It  is  worthwhile  noting  the  effects  that  various  groupings  can 
have  on  a  given  collection  of  fragments  and  how  this  can  influence  our 
mathematical  representation.  In  Table  I  we  present  a  sample  population 

Table  I.  SAMPLE  FRAGMENT  POPULATION  (Masses  in  Grams) 


1 

.106 

8 

.264 

15 

1.250 

22 

3.950 

2 

.110 

9 

.268 

16 

1.411 

23 

4.922 

3 

.115 

10 

.311 

17 

1.706 

24 

5.700 

4 

.123 

11 

.450 

18 

1.972 

25 

5.850 

5 

.151 

12 

.525 

19 

2.002 

26 

7.106 

6 

.172 

13 

.713 

20 

2.150 

27 

9.760 

7 

.195 

14 

.809 

21 

3.670 

28 

10.500 

22 


Figure  3.  Incomplete  Gaunma  Function. 


which  has  been  constructed  to  simulate  much  larger  populations 
typical  of  real  warhead  natural  fragmentation*  The  advantage  of 
using  such  a  small  population  is  that  it  is  easy  to  follow  in  detail 
various  possible  procedures.  Here  we  have  28  fragments  with  total 
mass  =  66.261g  and  average  mass  m  =  66.261/28=2 ,366g.  We  can 
represent  this  collection  graphically  by  drawing  28  vertical  lines  of 
height  unity  along  a  horizontal  mass  axis.  These  lines  will  be  more 
closely  spaced  for  smaller  m  values  and  become  farther  apart  as  m 
increases.  Alternatively  we  can  group  the  fragments  into  intervals 
centered  on  various  mass  values  and  count  the  number  in  each  group. 

If  we  choose  equal  size  intervals  each  one  mass  wide,  we  obtain 

Table  II.  As  we  see,  the  number  of  fragments,  N  ,  in  each  group. 


Table  II.  GROUPING  INTO  EQUAL  SIZE  MASS  INTERVALS 


Interval (g) 

0-1 

1-2 

2-3 

3-4 

4-5 

5-6 

6-7 

7-8 

8-9 

9-10 

10-11 

ne 

14 

4 

2 

2 

1 

2 

0 

1 

0 

1 

1 

nb(>) 

14 

10 

8 

6 

5 

3 

3 

2 

2 

1 

0 

nhC<) 

14 

18 

20 

22 

23 

25 

25 

26 

26 

27 

28 

fluctuates  erratically  for  the  heavier  groups  which  contain  only  a 
few  fragments.  This  is  typical  of  real  fragment  distributions  also. 

The  number  of  fragments  with  mass  greater  than  that  associated  with 
a  given  group  is  N^(:  ),  while  its  complement,  the  numbgr  with  pass 
less  than  that  of  higher  groups,  is  N^(<)  .  Of  courseN^O)  +  N"(<)  =  NT> 
the  total  number  of  fragments.  Both  the  cumulative,  N  O) ,  and  £ 
complementary  cumulative  numbers  appear  somewhat  smoother  than  N  ,  so 
we  might  expect  better  agreement  if  we  fit  a  smooth  mathematical 
function  to  cumulative  numbers  rather  than  to  the  frequency  N  , 

We  can  smooth  out  the  erratic  behavior  of  N  by  choosing  unequal  size 
intervals.  For  example,  the  choice  in  Table  III  gives  much  smoother 
behavior.  Of  course  other  choices  could  introduce  erratic  behavior 


Table  III.  GROUPING  INTO  UNEQUAL  SIZE  INTERVALS 


Interval (g)  0-.25 

.25-. 75 

.75-2 

2-4  4-6  6-10 

10-66.261 

NE  7 

6 

5 

4  3  2 

1 

NE(>)  21 

15 

10 

6  3  1 

0 

NEC<)  7 

13 

18 

22  25  27 

28 

24 


again.  For  example,  choosing  1.9S  instead  of  2.0  would  put  4  fragments 
in  the  third  group  and  5  in  the  fourth  group.  This  frequently  occurs 
in  practice  where  the  choice  of  mass  intervals  is  made  before  an 
experiment  is  carried  out,  perhaps  for  the  sake  of  uniform  reporting 
procedures.  Generally  speaking,  in  the  literature  the  mass  of  each 
fragment  is  not  reported  and  only  information  about  pre-chosen  groups 
is  given.  Consequently,  there  is  no  way  to  choose  new  groupings  which 
might  be  more  closely  represented  by  certain  mathematical  functions . 

Now  let  us  apply  the  frequency  distribution  dF=d(N/N  )  of  Equation 
(20)  with  y= m/p  and  B=y/p  where  p  is  the  density  of  the  metal  case. 

If  we^multiply  it  by  the  total  number  of  fragments,  N  ,  and  divide  by 
(l-e“  )  where  u  is  given  by  Equation  (27) ,  we  obtain 


dN 


N„ 


(l-e"U) 


exp 


dm  =  W(m)  dm  (33) 


for  the  number  of  fragments  in  the  infinitesimal  interval  from  m 

to  m+dm.  Note  that  Equation  (33)  becomes  Equation  (1)  above  for 

u-x”,  m  =0  and  £=3.  The  factor  (l-e~u)  insures  that 
o  v  J 


NT/(l-e'u 


exp 


NT/(l-e 


•u. 


/ 


u 


e  L  dt 


=  N„ 


(34) 


The  cumulative  number  of  fragments  is 


N(<) 


(35) 


while  the  complementary  cumulative  number  or  number  greater  than  m 
is  the  integral  from  m  to  in  (or  t  to  u)  . 


(36) 
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so  N(<)  +  N(>)  =  N  We  can  calculate  either  cumulative  number, 
however,  N(>)  is  the  number  used  to  evaluate  the  lethality  of  a 
warhead  or  the  vulnerability  of  a  target.  We  can  determine  N(>) 
experimentally  with  greater  accuracy  also,  since  the  number  and  size 
of  small  fragments  is  difficult  to  measure.  The  use  of  a  cutoff  mass, 
m  ,  of  sufficient  size  avoids  this  difficulty,  and  is  also  useful 
bicause  very  small  fragments  are  usually  not  lethal  anyway. 

In  practice,  we  do  not  deal  with  infinitesimals  such  as  dN  and  dm. 
Instead  we  use  fragment  groupings  such  that  N.  fragments  are  found  in 
the  itn  mass  group  of  width  dm4  centered  on  miss  m.. .  Generally 
speaking  also  in  fragmentationxwork  the  expected  n&mber  of  fragments 
is  large  and  approximately  equal  to  the  expected  number  of  breaks,  u. 

If  u  >  10,  then  l-e”u  1  and  may  be  neglected  in  Equations  (33) 
through  (36)  .  This  may  not  be  true  for  t  very  close  to  0  or  u  in 
Equations  (35)  and  (36),  but  the  use  of  finite  groupings  prevents 
this  from  happening.  This  is  why  u  can  be  taken  to  be  effectively 
infinite,  at  least  in  cases  of  natural  fragmentation.  For  a  narrow 
group  of  controlled  fragments  clustered  about  y  this  may  not  be  true 

for  small  £,  since  u  =  [(m  -m  )/y]1^  with  m  ,  m  and  y  all  of  about  the 
size.  For  finite  size  int&rvSls  Equation  (!$3)  8ecomes 


where 


Ni  =  nt  Wi 


(1-e 


(37) 

(38) 


for  the  itE  mass  group.  We  can  use  Equation  (38)  for  t.  instead  of 
t  in  Equations  (35)  and  (36)  to  take  account  of  finite  rather  than 
infinitesimal  size  mass  intervals. 


In  Equation  (37),  =  N../NT  is  the  calculated  probability  of 

finding  fragments  in  the  i  group.  The  probability  observed 

experimentally  is  W.E  =  N.E  /  N„,  where  N.E  is  the  number  found 
experimentally.  We1are  interested  in  obtaining  functional  represen¬ 
tations  of  experimental  data,  using  the  probability  formulas  we  have 
derived.  As  we  have  seen,  there  is  no  reason  to  require  that  £  be 
integer,  so  we  will  treat  it  as  an  adjustable  parameter.  In_addition, 
the  observed  average  mass  gives  us  only  a  rough  estimate  of  m  (and 
so  y  through  Equations  (30)  to  (32)).  We  will  use  this  as  an  initial 
guess  for  y  and  then  adjust  y  also  to  represent  the  data  better.  The 
cutoff  mass,  m  ,  will  determine  the  total  number  of  fragments,  N^,. 

For  example,  i8  Table  I  we  are  assuming  that  fragments  less  than  mQ  = 
0.1  g  are  of  no  interest  in  a  particular  application  because  they 
are  too  small  to  be  lethal.  If  m_  were  0.2  g  instead,  then  would 
be  21  instead  of  28.  If  m^  were  smaller  than  0.1  g  then  would 
be  larger  than  28.  Howeve?,  the  smaller  we  make  mQ,  the  mdre  uncertain 
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we  are  about  the  value  of  N  above  in  because  of  the  difficulty  of 
observing  very  small  fragments.  The°minimum  value  of  mQ  is  well 
above  molecular  mass  values,  since  we  are  considering  a  defect 
method  of  fragment  production.  Still  the  actual  number  of  dust-like 
fragments  emitted  after  a  detonation  is  undoubtedly  very  large 
and  experimentally  unknowable.  Since  this  number  is  also  of  no 
practical  interest  for  lethality  or  vulnerability,  we  will  adopt  the 
point  of  view  that  NT  is  fixed  by  a  choice  of  finite  mo>0.  The 
upper  mass,  m  ,  can  be  taken  to  be  the  total  mass,  NL, 'which  is 
usually  closeuto  the  unfragmented  case  mass .  For  natural  fragmentation 
this  will  make  u  effectively  infinite.  For  controlled  fragment  groups 
where  m  ,  p  and  m  may  be  close  to  each  other,  the  choice  may  be  more 
important  as  we  sHall  mention  later. 

If  our  only  interest  were  to  use  Equation  (36)  with  u  effectively 
infinite,  then  we  could  adopt  In  N.  (>)  as  our  model  function, 

1  1/ % 

considering  it  to  be  a  linear  function  of  the  variable  (nh-mo)  in 

In  N^)  =  In  NT  -  (p"A/A')  Cmi  -  mo)"'~  .  (39) 

For  fixed  4  we  could  adjust  p  and  so  the  slope.  For  fixed  p,  an 
adjustment  of  4  becomes  much  harder  since  we  are  changing  the 
independent  variable  as  well  as  the  slope.  Likewise,  simultaneous 
adjustment  of  4  and  p  is  too  laborious  to  consider.  The  logarithm  of 
Equation  (37)  is 


1/*, 


In  N.  =  In  [(NTAm.)/(£pi/A')]  +  (1-4)  In 


,  .1/4  ,  -1/4. 

-  ip  J 


(mi-mo)1/J'(40) 


which  is  not  a  linear  model  function,  except  for  4=1.  In  addition, 
if  u  is  not  effectively  infinite,  there  is  no  way  to  make  either  p  or 
4  appear  in  a  linear  fashion  in  a  model  function  related  to 
Equation  (36)  or  Equation  (37) . 


Here  we  will  use  the  least  squares  method  to  adjust  the  two 
parameters  p.  =  4  and  p~  =  p  which  appear  in  the  non-linear  model 
functions.  Equation  (36J  and  Equation  (37) .  The  function  to  be 
minimized  is  the  sum  of  the  squared  differences  between  the 
experimental  and  calculated  values  which  we  will  call  Sqd  (N  =  number 
of  groups) : 

N 


Sqd 


[Nt  (W.b  -  w.)]2 


(41) 


i=l 
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We  can  use  our  data  to  obtain  initial  guesses  for  the  parameters, 
namely,  p  =  Zq  and  p_0  =  u  .  If  these  guesses  are  reasonable,  we 
can  neglect  all  but  linear  correction  terms  in  an  expansion  of  W. 
about  the  point  (P1o,  p2o),  so  1 


(42) 


where  subscript  zero  means  evaluation  at  the  current  guess  point. 

Our  corrected  parameter  set  (Pk0+Ck  for  k=l,2)  becomes  the  new 
guess  point  in  an  iterative  calculation,  approaching  a  best  fit  as 
closely  as  we  please.  We  now  put  Equation  (42)  into  Equation  (41), 
set  the  derivatives  of  Sqd  with  respect  to  Z  and  p  equal  zero, 
divide  by  (-2N  ~)  and  obtain  the  normal  equations.  The  k““  equation 
is  1 


(43) 


where  k=l,2...  N  with  N  equal  to  the  number  of  adjustable  parameters 
(only  two  here) A  rearPangement  of  Equation  (43)  gives 


or,  in  matrix  notation 

QC  =  K  (45) 


where 
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(46) 


are  the  elements  of  the  symmetric  matrix  Q  of  dimension  N 


P 


K 


k 


x  N 


P 


and 


(47) 


are  the  elements  of  the  vector  K  with  the  correction  vector  C  =  Q  K. 
Qf  course  in  our  case  with  N  =  2,  the  solution  for  C  is  especially 
simple  since  Equation  (43)  consists  of  only  two  linear  simultaneous 
equations  for  the  unknowns  C.  and  C„,  More  elaborate  forms  of  the 
least  squares  method  could  also  be  dsed,  but  this  is  sufficient  for 
our  purpose  here. 


To  carry  out  this  procedure  we  need  derivatives  of  our  model 
functions  which  contain  factors  of  the  form  A  where  A  =  (m.-m )/y 
and  f  *  1/1  or  (—  -1).  We  recall  that  the  derivative  with  respect 
l  can  be  found  by  letting 


to 


Y  =  In  (Af)  =  f  In  A 


so 


and 


then 


dY  -f  d_ 

dl~  dH 


sr  ln  A  (af) 


d  W. 
1 

dt 


(y*)  j-i-a-t.) 

d  K.  r 

-  [Ki  > 


in  +  ue"u  in  u  /  (l-e*U)j 

J  -1  +  ti  +  ue"u  /  (l-e“u)J 


(48) 

(49) 

(50) 

(51) 

(52) 
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where  is  given  in  Equation  (37)  and  is  given  in  Equation  (38) . 

If  we  wish  to  optimize  a  fit  of  N. (>)  to  experimental  data,  the 
function  to  be  minimized  is 


N 

Sqd  =2|NiHc>)  "  Ni(>)  | 


(53) 


and  a  new  set  of  normal  equations  can  be  found  in  a  similar  way.  The 
required  derivatives  are 


A  computer  program  implementing  this  procedure  is  given  in  the 

Appendix.  For  cases  in  which  u  is  effectively  infinite,  all  terms 

involving  u  in  Equations  (51),  (52),  (54)  and  (55)  vanish  as  can 

easily  be  verified  by  the  use  of  L'  Hospital's  rule.  Provisions 

are  also  made  in  the  appendix  for  adjusting  either  i  or  p  alone, 

as  well  as  for  calculating  with  fixed  l  and  y,  using  either  finite 

m  and  m  or  with  m  =  0,  m  =  °°  as  in  Mott's  case, 
o  u  o  u 


C.  Applications 

Let  us  apply  our  procedure  to  the  sample  data  of  Table  I  as 
grouped  in  Tables  II  and  III.  For  example.  Table  IV  A  compares 
the  number  in  each  group  calculated  by  various  procedures  with  the 
number  found  experimentally,  N  .  In  this  Table  and  Table  V  m  was 
taken  to  be  zero.  If  we  took  m  =  0.1  g  in  Table  IV  we  could 
obtain  slightly  better  agreement  as  indicated  by  a  lowering  of  Sqd 
from  6.86  to  6.56,  accompanied  by  somewhat  changed  l  and  y.  However, 
this  is  good  enough  for  our  purpose,  which  is  to  show  that  the 
fitted  Sqd  is  much  lower  than  that  for  integer  £  values,  namely, 
33.53,  24.38  or  61.92.  For  the  fit,  m  was  taken  to  be  66.261  g 
so  that  u  =  (66 .261/ .95) -*  =  21.2  which  is  effectively  infinite, 

giving  N  approximately  equal  to  zero  for  the  last  group.  In  the 
Mott-type  calculations  m^  =  u  =  °°,  so  N  and  N(>)  are  both  exactly 
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zero.  This  is  true  because  the  exponential  factors  vanish  and 
dominate  other  factors  which  increase  without  limit.  Table  IV  B 
compares  calculated  values  with  experimental  values  of  N(>).  A  gain 
from  a  comparison  of  Sqd  values  we  see  the  value  of  adjusting  l  and  y. 
The  adjusted  Jl  =  2.03  is  almost  integer,  so  the  improvement  over 
l  =  2,  y  =  1.183  g  is  not  as  great  as  in  Part  A. 

Table  V  makes  the  same  comparisons  as  Table  IV  but  for  the 
sample  data  as  grouped  in  Table  III. 

Table  IV.  FRAGMENT  DISTRIBUTION  FOR  SAMPLE  IN  TABLE  II 


(Columns  marked  *1,  2  or  3  use  the  Mott  version.) 

A.  N  in  each  group  B.  N(>)  for  each  group 


Interval (g) 

ne 

N 

11=1 

1  =  2 

*<=> 

II 

04 

N(>) 

E  N  (>) 

*=> 

II 

H4 

II 

1  =  3 

0-1 

14 

13.51 

9.58 

9.50 

6.84 

14 

14.96 

22.67 

14.62 

9.48 

1-2 

4 

4.65 

6.28 

3.41 

2.04 

10 

9.54 

14.85 

9.08 

5.88 

2-3 

2 

2.18 

4.11 

1.90 

1.09 

8 

7.01 

9.73 

6.54 

4.40 

3-4 

2 

1.14 

2.70 

1.23 

.70 

6 

5.46 

6.38 

5.01 

3.53 

4-5 

1 

.64 

1.77 

.86 

.49 

5 

4.41 

4.18 

3.98 

2.95 

5-6 

2 

.38 

1.16 

.64 

.37 

3 

3.64 

2.74 

3.24 

2.52 

6-7 

0 

.23 

.76 

.48 

.29 

3 

3.05 

1.79 

2.69 

2.20 

7-8 

1 

.14 

.50 

.38 

.23 

2 

2.60 

1.18 

2.26 

1.94 

8-9 

0 

.09 

.33 

.30 

.19 

2 

2.23 

.77 

1.92 

1.73 

9-10 

1 

.06 

.21 

.25 

.16 

1 

1.94 

.51 

1.65 

1.56 

10+ 

1 

.00 

0 

0 

0 

0 

.14 

0 

0 

0 

l 

- 

1.39 

1 

2 

3 

- 

2.03 

1 

2 

3 

u(.g) 

- 

.95 

2.366 

1.183 

0.394 

- 

1.29 

2.366 

1.183 

0.394 

Sqd 

- 

6.86 

33.53 

24.38 

61.92 

- 

4.46 

106.43 

6.00 

62.00 
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Table  V.  FRAGMENT  DISTRIBUTION  FOR  SAMPLE  IN  TABLE  III 


(Columns  marked  £  =  1,  2  or  3  used  the  Mott  version.) 


A.  N  in  each  group  B.  N(>)  for  each  group 


Interval (g) 

NE 

N 

£  =  1 

£  =  2 

£  =  3 

n(>)e 

N  (>) 

£  =  1 

£  =  2 

£  =  3 

0-.25 

7 

6.24 

2.81 

6.58 

6.44 

21 

21.11 

26.56 

20.23 

14.16 

.25-. 75 

6 

5.15 

4.79 

4.75 

3.42 

15 

15.26 

22.67 

14.62 

9.48 

.75-2 

5 

5.44 

8.27 

4.67 

2.82 

10 

9.69 

15.66 

9.53 

6.14 

2-4 

4 

3.56 

6.66 

3.02 

1.71 

6 

5.47 

7.88 

5.70 

3.92 

4-6 

3 

1.68 

2.86 

1.47 

.85 

3 

3.22 

3.38 

3.58 

2.72 

6-10 

2 

1.44 

1.61 

1 .35 

.83 

1 

1.69 

.95 

2.08 

1.83 

10+ 

1 

.18 

0 

0 

0 

0 

.03 

0 

0 

0 

i l 

- 

1.72 

1 

2 

3 

- 

1.81 

1 

2 

3 

M(g) 

- 

1.41 

2 . 366 

1.183 

0.394 

- 

1.23 

2.366 

1.183 

0.394 

Sqd 

- 

4.40 

38.01 

6.55 

23.94 

- 

.98 

125.38 

2.56 

97.25 

As  we  see  the  fit  gives  better  agreement  than  the  use  of  integer  £  and 
appropriate  p.  If  we  ignore  £=1  (which  is  never  used),  only  one  case 
(for  £= 3)  gives  a  worse  fit  than  is  found  in  Table  IV.  We  expect 
this  since  the  data  in  Table  III  was  made  to  appear  smoothly 
decreasing  by  a  choice  of  mass  intervals.  For  £=1  this  choice 
makes  the  calculations  appear  very  erratic  since  a  maximum 
of  8.27  is  reached  for  the  third  mass  group.  The  probability  per 
unit  mass,  w(m)  in  Equation  (33),  decreases  steadily  as  m  increases 
for  any  ££l .  However,  W.  =  w.  Am.  in  Equation  (37)  is  the  product 
of  w.  with  Am.  so  that  ai*l  inci'easl  in  Am.  as  in  the  first  intervals 
of  T&ble  V  ca?l  more  than  offset  a  decrease  in  w. ,  resulting  in  a  net 
increase.  This  often  happens  to  experimental  d&ta  too,  when  we  use 
variable  size  groupings.  If  all  mass  intervals  are  equal,  then  the 
calculated  numbers  will  always  decrease  as  m  increases.  However, 
experimentally,  erratic  behavior  can  occur  even  for  equal  size 
intervals  because  of  poor  statistics.  In  Table  V  A  with  £=1, 

N.  =  (28/2.366)  Ann  exp  (-m./2.366).  This  gives  N  =  2.81  for  Ann  = 
.25,  m  =  .125,  etc.  as  sho&n.  For  larger  £  values  this  difficulty 
tends  to  disappear. 
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Now  let  us  examine  some  experimental  data  obtained  from  a  real 
munition.  In  the  Netherlands  Lindeijer  and  Liemans18  published  data 
in  the  open  literature  for  a  105  mm  shell  which  is  presented  here  as 
the  experimental  values  in  Table  VI.  An  unusual  feature  of  their 
experiment  was  the  large  number  of  tiny  fragments  they  collected  in 
the  sub- gram  category  compared  to  the  816  individually  weighed  and 
counted  (and  recorded  in  Table  VI).  They  ve.re  able  to  collect  2,195 
tiny  fragments  weighing  215  grams  altogether,  compared  to  the 
816  recorded  here  which  totalled  457  grams.  The  total  number  of 
fragments  in  Table  VI  is  N„  =  2,228  weighing  a  total  of  NL.  =  8,561  g, 
giving  m_v  =  3.84  g.  This  population  is  typical  of  that  reported 
by  other  groups  in  similar  experiments.  Since  the  original  case 
mass  weighed  about  9,000  grams,  about  2.5  per  cent  was  lost, 
presumably  as  a  very  large  number  of  very  small  fragments  lost 
despite  the  extraordinary  efforts  of  these  authors.  Since  the  actual 
number  of  small  fragments  emitted  (or  created  by  secondary  breakup 
during  recovery)  cannot  be  determined  (although  it  must  be  extremely 
large),  we  see  the  need  for  invoking  a  fixed  value  of  N^,  determined 
by  a  non- zero  value  of  mQ. 

In  Table  VI  we  have  divided  the  Sqd  values  by  the  total  number  of 
fragments  to  keep  these  values  reasonably  small.  From  a  comparison  of 
these  values  it  is  clear  that  using  "best"  %  and  m  can  give  a  better 
overall  representation  of  experiment  than  the  use  of  integer  l  and 
prescribed  u  =  m.v/(£!).  To  avoid  crowding  in  Table  VI  we  have 
omitted  columns  Tor  values  found  using  £  =1  (worse  anyway)  and  have 
rounded  calculated  fragment  numbers  to  integer  values,  at  least 
when  they  exceed  unity.  The  Sqd  values  could  have  been  improved  a 
few  more  percent  in  the  "best"  calculations  by  using  mQ>0,  but  this  is 
not  shown  here  for  the  sake  of  uniformity  in  comparisons.  The  low 
values  found  for  the  heavest  groups  reflect  the  infinite  (or  effectively 
infinite)  upper  mass  limits  used.  Although  this  is  required  in  the 
Mott-Type  calculation,  improvements  in  the  "best"  calculation  can 
also  be  made  by  using  the  mass  of  the  heaviest  fragment  found  instead 
of  the  case  mass.  Note  that  the  "best"  l  values  are  not  very  different 
from  m  / (2 ! )  =  1.92  g,  although  the  "best"  SL  values  are  notably 
lower  than  i=2,  giving  better  agreement  for  the  more  populated  groups 
and  worse  for  those  with  only  a  few  fragments.  When  we  estimate 
effectiveness  for  cases  where  fragments  less  than  ten  grams  are 
important,  it  is  presumably  more  important  to  have  a  better  representa¬ 
tion  of  the  more  populated  groups,  since  there  are  so  few  heavy 
fragments.  If  these  lighter  fragments  are  not  important  for  a  particular 
application,  we  can  always  increase  m  ,  reduce  N  ,  adjust  Si  and  m 
and  again  obtain  "best"  values  for  th$  remainingpopulation.  Similar 
calculations  have  been  carried  out  for  other  real  munitions  (both 
larger  and  smaller  caliber)  and  similar  results  have  been  obtained. 


18  E.  W.  Lindeijer  and  J.  S.  Leemans,  Explosivstoffe  16_  (7) ,  145 
(1968). 
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Instead  of  using  any  type  of  analytical  formula,  one  can 
of  course  use  the  experimental  data,  especially  in  this  age  of  the 
digital  computer.  However,  for  some  purposes,  analytical  formulas 
are  desirable.  In  addition,  the  smoothing  which  results  from 
the  use  of  a  fitted  analytical  formula  should  better  represent  the 
average  values  and  average  effectiveness  of  a  large  number  of  munitions, 
the  case  we  are  trying  to  calculate.  Use  of  experimental  data  with 
poor  statistics  giving  erratic  behavior  in  heavier  mass  groups  (and 
certain  angular  zones)  could  give  poorer  effectiveness  estimates  than 
use  of  smooth  functions.  Presumably  the  data  would  become  smooth 
if  only  we  could  afford  to  do  enough  experiments.  Sometimes  a  few 
larger  fragments  might  be  neglected  anyway  because  they  are  too  slow 
to  be  of  importance. 

Figure  4  compares  the  curves  for  N(>)  from  Table  VI.  The  "best" 
fit  curve  is  virtually  indistinguishable  from  the  experimental  curve 
on  the  scale  of  this  graph. 


Table  VI.  FRAGMENT  DISTRIBUTION  FOR  THE  NATURAL  FRAGMENTATION 


Interval (g) 

A. 

ne 

N  in  - 

N 

OF  A  105  MM 

each  group 

a  =  2  a  =  3 

SHELL 

B 

nc>)e 

.  N(>) 

N  (>) 

for  each 

1=2 

group 

l  =  3 

0-1 

816 

804 

683 

544 

1412 

1421 

1337 

895 

1-2 

464 

414 

271 

175 

948 

938 

921 

598 

2-3 

234 

257 

162 

97 

714 

689 

712 

469 

3-5 

270 

282 

190 

109 

444 

471 

526 

360 

5-10 

233 

220 

203 

118 

211 

233 

309 

235 

10-15 

97 

53 

89 

55 

114 

104 

174 

155 

15-20 

49 

15 

47 

32 

65 

53 

109 

113 

20-25 

28 

5 

28 

21 

37 

29 

73 

87 

25-30 

13 

2 

17 

15 

24 

17 

51 

70 

30-35 

5 

.51 

12 

11 

19 

10 

36 

57 

35-40 

1 

.18 

8 

8 

18 

6 

27 

48 

40-45 

4 

.07 

6 

6 

14 

4 

20 

40 

45-50 

4 

.03 

4 

5 

10 

3 

15 

35 

50+ 

10 

.00 

0 

0 

0 

.00 

0 

0 

l 

- 

1.27 

2 

3 

- 

1.68 

2 

3 

y 

- 

2.18 

1.92 

.64 

- 

1.91 

1.922 

.64 

Sqd/NT 

- 

3.33 

30.38 

97.67 

- 

1.25 

13.69  ; 

210.60 
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Finally,  let  us  consider  an  application  to  controlled 
fragmentation.  It  is  obvious  that  for  particular  targets  or  target 
classes  there  will  be  an  optimum  fragment  size  (also  dependent 
somewhat  on  engagement  conditions) .  Smaller  fragments  are  wasted 
since  they  cannot  sufficiently  damage  the  target,  while  fragments 
which  are  too  large  will  kill  the  target  but  inefficiently  since 
they  will  of  necessity  be  fewer  in  number,  so  less  likely  to  hit 
without  special  aiming  devices.  The  ultimate  in  control  is  the 
preformed  fragment  which  is  often  used  when  launch  stresses  are 
relatively  low  as  for  missile  warheads.  However,  gun-launched 
warheads  require  strong  casings  which  do  not  permit  weakening  the 
shell  by  grooving.  For  such  shells  other  techniques  have  been 
devised  which  permit  undiminished  case  integrity  coupled  with  some 
degree  of  fragment  size  control.  Examples  are  the  use  of  electron 
beam  scoring,  selective  carburization  of  the  steel  and  grooving  of 
the  explosive  at  the  shell  wall  (with  or  without  a  liner) .  At  BRL 
Meissner  and  Kineke19  have  compared  various  techniques  with  each 
other  and  with  natural  fragmentation.  For  example,  they  used  a 
cylindrical  steel  case  10.16  cm  long  with  outer  diameter  7.6  cm  and 
wall  thickness  0.24  cm,  filled  with  Composition  B  explosive  to  make 
such  comparisons.  Table  VII  compares  the  natural  fragmentation  they 
observed  with  values  calculated  by  various  methods .  Clearly  a  better 
representation  can  be  obtained  by  using  best  values  of  l  and  y. 

Since  the  mass  of  the  largest  fragment  was  recorded  and  is  known  to 
be  less  than  14  g,  this  value  was  used  for  m^. 

In  other  experiments  Meissner  and  Kineke  have  grooved  the  case  or 
the  explosive  in  an  effort  to  control  fragment  size.  In  particular 
they  used  22  longitudinal  grooves  evenly  spaced  together  with  2 
circumferential  grooves,  each  one  1.27  cm  from  either  end.  For  the 
sake  of  illustration  we  will  discuss  the  results  they  obtained  using 
a  grooved  explosive.  As  we  see  from  Figure  5  and  Table  VIII,  they 
succeeded  in  producing  41  fragments  with  mass  greater  than  13  g, 
compared  to  only  one  such  fragment  when  no  effort  was  made  to  control 
fragment  size  (Table  VII).  These  rod-like  fragments  centered  about 
a  mass  near  17  g.  In  addition,  they  produced  about  21  end  fragments 
in  a  group  which  centered  near  6  or  7  g.  The  rest  of  the  case  mass 
seems  to  be  in  a  naturally  fragmented  group  with  most  of  the  mass 
below  5  g.  Column  one  of  Table  VIII  (N  )  gives  the  number  observed  in 
each  mass  interval.  This  column  has  been  divided  into  the  thrge  groups 
in  the  next  three  columns,  using  Figure  5  $s  a  gui^e.  Here  N,  is  the 
natural  part  of  the  fragmentation  while  N9E  and  are  the  partially 


19  R.  Meissner  and  J.  Kineke,  report  to  be  published. 
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MASS  (g) 


Z 


v 


Table  VII.  FRAGMENT  DISTRIBUTION  FOR  THE  NATURAL  FRAGMENTATION 

OF  A  CYLINDER 

A.  N  in  each  group  B.  N(>)  for  each  group 


Interval (g) 

ne 

N 

CM 

11 

o* 

A  =  3 

n(>)e 

N(>) 

CM 

II 

o* 

a  =  3 

0-1 

216 

205 

156 

110 

235 

252 

226 

146 

1-2 

65 

81 

54 

32 

170 

154 

136 

89 

2-3 

47 

42 

30 

17 

123 

107 

96 

66 

3-4 

44 

24 

19 

11 

79 

78 

73 

52 

4-5 

20 

14 

13 

8 

59 

58 

57 

43 

5-6 

17 

9 

10 

6 

42 

44 

46 

37 

6-7 

21 

6 

7 

4 

21 

34 

37 

32 

7-8 

5 

4 

6 

3 

16 

25 

31 

28 

8-9 

8 

3 

4 

3 

8 

19 

26 

25 

9-10 

3 

2 

4 

2 

5 

14 

22 

22 

10-11 

0 

1 

3 

2 

5 

10 

19 

20 

11-12 

3 

.83 

2 

1.7 

2 

6 

16 

18 

12-13 

1 

.59 

2 

1.5 

1 

3 

14 

17 

13-ra 

u 

1 

.42 

0 

0 

0 

1 

0 

0 

% 

- 

1.35 

2 

3 

- 

1.83 

2 

3 

y 

- 

1.24 

1.05 

.35 

- 

1.46 

1.05 

.35 

Sqd/NT 

- 

2.61 

10.97 

33.21 

- 

2.91 

8.17 

45.06 

Note:  mu  = 

14g 

for  best  calculations 

m  = 
u 

infinity 

for  l 

=  2  or 

controlled  parts.  Since  the  controlled  parts  exhibit  maxima,  we 
must  use  JI<1  in  Equation  (37)  to  represent  them.  This  may  be  seen 
by  examining  the  mass -dependent  part  of  Equation  (37) , 
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where  has  been  included  to  make  it  equal  to  w  in  Equation  (33)  * 

For  £>¥,  w  decreases  smoothly  from  infinity  at  m  =  m  to  zero  as  m 
becomes  large,  as  we  have  seen  in  our  representation?  of  natural 
fragmentation.  For  Jt=l  we  have  a  simple  exponential,  finite  at 
m  =  m  .  For  0<£<1,  w  vanishes  for  m  =  m  or  infinity,  but  has  a 
maximum  in  between.  0 

In  Table  VIII  the  columns  marked  N^,  ^  and  give  the  fragment 


TABLE  VIII.  PARTIALLY  CONTROLLED  FRAGMENTATION  OF  A  CYLINDER 

(blanks  mean  zero) 


Interval 

ne 

N1E 

»2e 

"sB 

N 

N1 

N2 

N3 

n(>)e 

N  (>) 

o-i (g) 

63 

63 

61.22 

61.22 

94 

89.04 

1-2 

15 

15 

18.33 

18.33 

79 

70.71 

2-3 

13 

13 

5.78 

5.78 

66 

64.93 

3-4 

3 

3 

1.87 

1.87 

63 

63.06 

4-5 

6 

1 

5 

4.32 

.61 

3.71 

57 

58.74 

5-6 

5 

5 

7.00 

7.00 

52 

51.79 

6-7 

8 

8 

5.86 

5.86 

44 

45.88 

7-8 

3 

3 

3.20 

3.20 

41 

42.68 

8-9 

1.24 

1.24 

41 

41.44 

9-10 

.36 

.36 

41 

41.08 

10-11 

41 

41.08 

11-12 

41 

41.08 

12-13 

41 

41.08 

13-14 

1 

1 

1.97 

1.97 

40 

39.11 

14-15 

6 

6 

4.61 

4.61 

34 

34.50 

15-16 

4 

4 

6.09 

6.09 

30 

28.41 

16-17 

10 

10 

6.48 

6.48 

20 

21.93 

17-18 

3 

3 

6.03 

6.03 

17 

15.90 

18-19 

5 

5 

5.06 

5.06 

12 

10.84 

19-20 

4 

4 

3.89 

3.89 

8 

6.95 

20-21 

5 

5 

2.77 

2.77 

3 

4.18 

21-22 

1 

1 

1.83 

1.83 

2 

2.35 

22-23 

0 

0 

1.14 

1.14 

2 

1.21 

23-24 

1 

1 

.66 

.66 

1 

.55 

24-25 

0 

0 

.36 

.36 

1 

.19 

25-26 

1 

1 

.19 

.19 

0 

0 
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Figure  5.  Partially  Controlled  Fragmentation  of  a  Cylinder. 


numbers  calculated  using  adjusted  £  and  y  values*  For  N  the  total 
number  of  fragments  was  taken  to  be  95  (as  in  N  E)  while  m  =0,  giving 
£=^.04  y=.82  g  and  Sqd/N  =.77.  For  N2  the  total  number  N°  =21  (as  in 
N2  ).  In  addition,  the  values  m  =4  g  and  m  =10  g  were  usea,  giving 
£-.53,  y=2.47  g  and  Sqd/T^j=.56.°  Here  the  ?um  of  mQ  and  y  gives  an 
indication  of  the  peak  position.  Using  mn>0  seems  to  give  a  better 
representation  than  m  =0  here.  For  N_,  NT_=41  and  m  =13  g,  giving 
£=.54,  y=5.18  and  Sqd/N  =.89.  The  calculated  value?  for  these 
groups  are  shown  by  the  solid  lines  in  Figure  5.  The  agreement 
appears  to  be  reasonable.  The  last  two  columns  in  Table  VIII  are  the 
complementary  cumulative  distributions  obtained  by  subtraction, 
starting  with  the  total  number  observed,  N  =157,  or  with  the  total 
number  calculated  N,j=150.26.  Again  the  agreement  is  reasonably  good. 


Strictly  speaking,  we  should  carry  out  a  simultaneous  fit  of 
three  theoretical  curves  to  the  data,  using  m  and  N  as  well  as  £ 
and  y  for  each  curve  as  adjustable  parameters?  This  more  complicated 
procedure  (or  one  similar  to  it)  is  actually  used  in  certain  types 
of  spectroscopy  where  the  statistics  are  extremely  good  in  order  to 
extract  information  about  partially  resolved  lines.  However,  in  the 
present  case  where  the  statistics  are  so  poor,  the  procedure  used 
seems  adequate. 


III.  SUMMARY 

We  have  seen  that  Mott's  distribution  is  a  special  case  of  more 
general  probability  distributions  which  are  widely  used  to  describe 
phenomena  in  various  fields .  In  spite  of  this  it  can  be  given  a 
rational  basis  in  terms  suited  to  fragmentation,  namely,  the 
random  distribution  of  defects  in  solids.  Methods  of  grouping 
fragment  populations  have  been  discussed,  and  a  general  probability 
distribution  has  been  applied  to  several  collections  of  experimental 
data.  In  particular,  it  has  been  shown  that  this  approach  can  be 
used  not  only  for  natural  fragmentation,  but  also  for  controlled 
fragmentation. 
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APPENDIX 


This  computer  code  for  the  Hewlett-Packard  9845B  offers  a  number 
of  options  in  fitting  Equations  (36)  and  (37)  of  the  text  to 
experimental  fragment  size  data.  These  should  be  self-evident  from 
the  various  questions  asked  by  the  interactive  display  statements 
which  occur  throughout  the  program.  Since  the  first  part  of  the 
program  aims  to  optimize  a  fit  to  the  number  in  each  group,  the 
output  subroutine  (lines  1470  to  1730)  always  prints  N  values  but 
permits  N(>)  values  to  be  printed  if  desired,  even  though  these  have 
not  been  optimized.  Similarly,  in  the  output  subroutine  for  part 
two  (lines  2970  to  3220),  N(>)  is  always  printed,  but  N  is  an  option. 
This  is  particularly  useful  when  the  program  is  used  without  adjusting 
the  parameters.  This  occurs  in  two  sections,  one  starting  on  line 
1030  in  which  the  values  of  m  and  m  specified  at  input  are  used, 
the  other  starting  on  line  24^0  in  wKich  m  =0  and  m  =°°. 

O  U 
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10  REM  FRAG.  FIT  TWO- PARAMETER  FUNCTIONS  TO  DATA.  FIRST  PART  FITS  NUMBER  VS 
MASS  GROUP.  SECOND  PART  FITS  COMPLEMENTARY  CUMULATIVE  NUMBER . 

20  REM  INPUT:  MM  I  N  =  CUTOFF  MASS , NTOT=TOTAL  NUMBER  OF  FRAGMENTS ,  NG=NUMBER  OF  MA 
GROUPS  ,  MIJP  C  I  >  =UPPER  LIMIT  OF  I  TH  GROUP ,  N  <  I  >  =NUMBER  I  N  GROUP  I . 

30  REM  P 1  =  L  AND  F'2=MU  ARE  ADJUSTABLE  PARAMETERS. 

40  OPTION  BASE  1 

59  PRINTER  IS  0 

60  11=0 

7  O  D I M  M  <  5  0  > ,  W  <  5  0  >  ,  D  <  5  0  > ,  D  e  v  <  5  0  )  ,  M  u  p  <  5  0  > ,  N  <  5  0  >  ,  D  rn  <  50>,  N  g  m  c  <  5  0  >  ,  N  g  m  e <  5  0 ) ,  V  < 

£  0  DIM  D  v  <50  > 

SO  IMAGE  2D, 5<8D. 2D) 

100  DISP  "MM  IN, NTOT, NG" ; 

110  I N  P  U  T  M  rn  i  n ,  H  t  o  t  ,  N  g 

120  PRINT  "  MM  I  N=  "  ;  Mm  i  ri ;  "NTOT="  ;  Ntot ;  "  N  G  =  "  ;  N  g 

1  30  PR  I  NT  "  I  MUP < I >  N  < I )  NGME < I > " 

1 4 O  F OR  1  =  1  T 0  N g 

150  DISP  " MUP <  I  ,  N <  I  >  "  ; 

160  INPUT  Mup < I > , N < I > 

170  IF  I>1  THEN  230 
1  8  O  D  rn  <  1  >  =  M  u  p  <  1  >  -  M  rn  i  n 

1  9  O  M 1  1  >  =  <  M  u  p  <  1  >  +  M  rn  i  n )  ■  2 

200  D< 1 >=N< 1 }  Ntot 

210  N  g  rn  e  <  1  >  =  N  t  o  t  -  N  <  1  ) 

220  GOTO  270 

2  3  O  D  rn  <  I  )  =  M  u  p  <  I  >  -  M  u  p  <  I  -  1  > 

2  1 G  M  <  I  >  =  <  M  u  p  <  I  >  +  M  u  p  <  I  -  1  )  >  /  2 

250  D  <  1  >=N<  I  >-'Nt  ot 

2  6  0  N  g  rn  e  <  I  )  =  N  g  rn  e  <  I  - 1  >  -  N  <  I  > 

2  70  PR  I  NT  US  I NG  90 ; I , Mup  < I > , N < I > , Ngme < I > 

280  NEXT  I 

296  DISP  "FIT  FRAGMENT  NUMBER  IN  EACH  GROUP?  <0  IF  NO,  1  IF  YES:-": 

300  INPUT  Ind 

318  IF  I nd=0  THEN  GOTO  1740 

320’  PRINT  "FIT  EXPRESSION  FOR  NUMBER  IN  ITH  GROUP  TO  DATA" 

330  PRINT 

340  DISP  "Pl=?  P 2  =  ?  0,0  MEANS  SKIP  THIS  PART  "; 

350  INPUT  PI , P 2 

360  IF  P 1 =0  THEN  GOTO  730 

370  PR  I  NT  "ADJUST  P 1  =  " ; P 1 ; "  AND  P2=" ; P2 

380  Q 1 1 =  Q 1 2  =  Q  2  2  =  K 1 =  K  2  =  S  q  d  =  0 

390  FOR  1=1  TO  Ng 

4  0  O  A  =  <  M  <  I  >  -  M  rn  i  n  >  /  P  2 

410  X  =  A ■■'■•<  1  '•PI  > 

4  2  O  B  =  <  M  u  p  <  N  g  >  -  M  m  i  n  >  P  2 

430  I J  =  B '  < 1  PI  > 

440  T 1 =  E  X  P  < - X ) 

450  T  2  =  E  X  P  < - U  > 

460  T  3  =  U  *  T  2  •'  <  1  -  T  2  > 

470  W  <  I  >  =Drn  <  I  )  *X*T  1  ✓  <  P 1  *F'2*A*  <  1 -T2  )  '> 

480  Nqrnc  <  I  >  =  N t  ot  *  <  T  1  - T 2  >  /  <  1  -T2  > 

490  D 1 =  W < I ) *  < - 1 -LOG  < X  >  *  < 1 -X  >  +T3*L0G <  U > )  -  P 1 

500  D 2  =  W  <  I  >  *  <  -  1  +X  +  T3  >  <  P  1  *P2  > 

510  01 1 =Q 1 1+D1" 2 

520  Q 1 2=0 1 2+D 1 *B2 

530  Q22  =  Q22  +  D2-'-2 

540  Dev < I )=D< I >-W< I > 

5  5  0  K 1  =  K 1  +  D 1  *  D  e  v  <  I ) 

560  K2=K2+D2*Dev< I ) 

570  Sqd=Sqd+ < Dev < I > *Nt ot > A2 
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580  NEXT  I 

590  PRINT  "SQDB" 5  Sqd 

600  Den=Ql 1*Q22-Q12A2 

610  Cl=<Kl*Q22-K2*Q12>2D*n 

620  C2*<K2*Q1 l-Kt*Q12)/D*n 

630  Rm=MAX<ABS<Cl/Pl ) ,  ABS<C2-'P2>  ) 

640  IF  Rm<*. 2  THEN  670 
650  Cl=.2/Rm*Cl 
660  C2  = .  2  ARm*C2 

679  IF  (RBStCl  )<  .  001  )  AND  <ABS<C2X  .  001 >  THEN  720 

680  P1=P1+C1 

690  P2=P2+C2 

700  print  *’ p i ■=  ";Pi;"P2*  ";P2 

710  GOTO  380 

720  GGSUB  Output  1 

730  BISP  "Pl=?  P2=?  ADJUST  PI  ONLY.  0,0  TO  SKIP"; 

740  INPUT  PI, M2 

750  IF  P1=0  THEN  GOTO  1030 

760  PRINT 

770  PRINT  WITH  P2=";M2;",  ONLY  ADJUST  P1  =  ";P1 

780  011 =  K 1 =Sqd=0 

799  FOR  1=1  TO  Ng 

800  A=. (M<  I  >-Mmi  nX'M2 
810  X=A ' < 1 /P 1 > 

828  U=  <  «■!  Mup  (  Ng  >  -Mm  i  n  >  /M2  >  A  <  1  -'PI  "> 

830  T2=EXPC-U) 

840  T3=U*T2/< 1-T2) 

850  W<  I  >=Dm<  I  )*X*EXP<-X>/'<P1*M2*A*<  1-T2) ) 

860  Ngmc  <  I  >=Ntot*<EXPC-X)-T2)-''<  1-T2) 

870  It  1  =W  <  I )  * <  - 1  -LOG  < X )  *  <  1  -X  )  +  T3*L0G  <  U )  )/Pl 

880  01 1=01 1+D1A2 

890  Dev< I )=D< I >-W< I > 

9 6 0  Sqd=Sqd+<Dev<  I  )*Ntot  ) A 2 

9 1 @  K 1 =K 1 +B 1 *D*m  < I > 

920  NEXT  I 
930  PRINT  " SQD= " ; Sqd 
940  C1=K1-'Q11 

950  Rm=ABSCCl-'Pl  > 

960  IF  R m <  =  . 2  THEN  980 

970  C 1  = .  2  R  m  *  G 1 

980  IF  ABSCC1X.001  THEN  1020 

990  P 1  =  P 1  +  C 1 

1000  PRINT  "P1=H;P1 

1010  GOTO  780 

1020  GOSUB  Output  1 

1030  IlISP  "Pl=?  P2  =  ?  BOTH  KEPT  CONSTANT.  0,0  TO  SKIP" 

1040  INPUT  Ml,  M2 

1050  IF  M 1 =0  THEN  GOTO  1170 

i860  S  q  d  =  0 

1070  PRINT  "  USE  CONSTANT  P1=";M1;"  AND  P2=";M2 

1080  FOR  1=1  TO  Ng 

1090  X=  <  <  M <  I  )  -Mm  i  n  >  -'  M2  >  X  1  AM 1  > 

1  1 00  U=  X  Mup  < Ng )  -Mm  i  n  > ^M2 >  X  1  /M 1  ) 

1110  T2  =  EXP  < -U ) 

1 1 20  W < 1 >  =Dm < I >  *X*EXP < -X  >  X  M 1  *  <  M< I > -Mm i n ) *  <  1 -T2 >  > 
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1130  Ngmc < I >=Ntot*<EXP<-X>-T2>/( 1-T2) 

1140  DeuCI )*B< I )-W<I > 

1150  NEXT  I 

1160  GOSUB  Output  1 

1170  DISP  "PI*?  P2»?  ADJUST  P2  ONLY.  0,0  TO  SKIP"; 

1180  INPUT  M1,P2 
1190  IF  M 1 =0  THEN  1450 

1200  PRINT  "  WITH  PI*"; Ml; " , ONLY  ADJUST  P2*";P2 

1210  Q22=K2*Sqd=0 

1220  FOR  1=1  TO  Ng 

1230  X*  < < M  < I > -Mm i n  > 'P2  > < 1 '  M 1 > 

1240  U*  <  (. Mup (. Ng  > -  Mm  i  n >  /P2 )  A (  1  /M 1 ) 

1250  T2  =  EXP<  -U ') 

1260  T3=U*T2/<1-T2> 

1270  W< I ) =Dm ( I >*X*EXP<-X>/<M1*<M< I >-Mmin>*< 1-T2>  > 

1280  Ngmc  < I ) =Nt ot * < EXP < -X ) -T2 ) ✓ < 1 -T2 ) 

1290  D2=N< I >*<-l+X  +  T3>-'<Ml*P2> 

1300  022  =  022  +  02''"  2 

1310  Deu< I >=D< I )-W(I ) 

1320  Sqd=Sqd+<Dev<I)*Nt ot ) 2 
1330  K2=K2+D2*Dew< I > 

1340  NEXT  I 
1350  PRINT  " SQD* " ; Sqd 
1360  C2=K2/Q22 

1370  Rm=AES<C2/P2) 

1380  IF  Rm<  = . 2  THEN  1400 

1390  C2  = .  2--'Rm*C2 

1400  IF  RBS  <  C2  ')  <  .  00 1  THEN  1440 

1410  P2=P2+C2 

1420  PRINT  " P2* 11 ;  P2 

1430  GOTO  1210 

1440  GOSUB  Output  1 

1450  GOTO  1740 

1460  REM  END  OF  PART  ONE 

1470  Output l:  Sqd=0 

1480  PRINT 

1490  PRINT  "  **NUMBER  IN  ITH  GROUP**  " 

1500  PRINT  "I  M  <  I  >  EXPER  CALC  DEVIATION" 

1510  FOR  1=1  TO  Ng 
1520  Tw  =  Ntot*W< I ) 

1530  Td=Ntot*D<I> 

1540  Tdeo=Nt ot *Deu ( I  > 

1550  Sqd=Sqd+Tdeu^2 

1560  PRINT  USING  90 ; I , M < I > , Td , Tw , Tdev 

1570  NEXT  I 

1580  PRINT  "SQD=" ; Sqd 

1590  DISP  "PRINT  NUMBER  WITH  MASS  GREATER  THAN  M  ?  <0,1>"; 

1600  INPUT  Ind 

1610  IF  lnd=0  THEN  GOTO  1720 

1620  Sqd=0 

1630  PRINT 

1640  PRINT  "  **NUMBER  WITH  MASS  GREATER  THAN  M**" 

1650  PRINT  "I  M  < I >  EXPER  CALC  DEVIATION" 
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1660  FOR  1*1  TO  Ng 
1670  Dev=Ngmc  <  I)-Ngiwe<  I  > 

1680  Sqd  =  S.qd+DevA2 

1690  PRINT  USING  90; I , M< I > , Ngme< I > , Ngmc < I >  ,  Dev 

1700  NEXT  I 

1710  PRINT  " SQD* " ; Sqd 

1720  PRINT 

1730  RETURN 

1740  PRINT 

1750  PRINT  "FIT  EXPRESSION  FOR  NUMBER  WITH  MASS  GREATER  THAN  M<.  I  .'  Tu  HATH 
1760  PRINT 

1770  .DISP  "Pl=?  P2  =  ?  <0,0  MEANS  SKIP  THIS  PART > " ; 

1780  INPUT  P 1 , P2 

1790  IF  P 1 *0  THEN  GOTO  2160 

1800  PRINT  "  ADJUSTING  PI*";  Pi;  "  AND  P2=";P2 

1810  Q 1 l=Q12=Q22=Kl*K2*Sqd*0 

1820  FOR  1=1  TO  Ng 

1830  A* <M< I )-Mm i n)/P2 

1840  X  =  A"< l'Pl  >  ' 

1850  U=  <  <  Mup  <  Ng  )  -Mm  i  n  >  -"P2  > A  <  1  P 1  > 

I860  T1=EXP<-X) 

1870  T2  =  EXP  < -U ) 

1830  T3=  <  1  -T 1  >  •/  <  1 -T2  > 

1890  W< I >  =  Dm< I >*X*T1/<P1*P^*A*< 1-T2> > 

1  900  Ngmc  <  I  )  =Nt  ot-  *  <  T 1 -T2  >  <  1 -T2  > 

1910  Dl=Ntot*<X*Tl*L0GOO-T3*U*T2#LOG(U))'Pl-'<l-T2J 
1920  D  2  =  N  t  o  t  *  <  X  *  T 1  -  T  3  *  U  *  T  2  >  P  1  -•  P  2  '  <  1  -  T  2  > 

1930  Q11=Q11+D1A2 

1940  Q  12  =  0 1 2  +  D 1 *D2 

1950  Q22=Q22+D2A2 

1 960  Dev < I ) =Ngme  < I ) -Ngmc < I > 

1970  Sqd=Sqd+Dev< I > A2 
1980  K 1 =K 1 +D 1 *Dev < I > 

1990  K2=K2+D2*Dev ( I ) 

2000  NEXT  I 

2010  PRINT  " SQD= " ; Sqd 

2020  Den=Ql 1*022-012 A2 

20  30  C  1  =  K  1  *Q22-K2*0 1 2 >  -  Den 

2040  C2=<K2*Q1 1-K1*Q12> -  Den 

2050  REM  LIMIT  SIZE  OF  CORRECTIONS 


2  0  6  0 

Rm=MAX<ABS<Cl 'PI > , ABS<C2  'P2>  > 

2070 

IF  Rm<  = . 2  THEN  2100 

2  0  3  0 

Cl  =  . 2 'Rm*C  1 

2  0  9  0 

C  2  = .  2  •'  R  m  *  C  2 

2  1  0  0 

IF  <ABS<C1 X . 001 >  AND  <ABS<C2> 

<  .001 

>  THEN  2150 

2110 

P1=P!+C1 

2120 

P2=P2+C2 

2130 

print  "Pi=  " ; P l ; " P2=  ";P2 

2140 

GOTO  1810 

2150 

GO SUB  Out  put 2 

2160 

HISP  "PI*?  P2=?  ADJUST  PI  ONLY 

.  0,0 

MEANS  SKIP"; 

2  1  70 

INPUT  PI,  M2 

2  i  36 

IF  P 1 =0  THEN  GOTO  2460 

2  1  9  0 

PRINT  "  ADJUST  Pi="; 

Pi;" 

ONLY  WITH  P2=";M2 

2  2  0  0 

Q 1 1  =K  1 =Sqd=0 

2210 

FOR  1=1  TO  Ng 

2  2  2  0 

A=<M< I >-Mmi n)'M2 

2  2  3  0 

X  =  H-<  1.--P1  ) 
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2240  U  =  ( < Mup(Ng) -Mm i  n)^M2 1  'PI  ) 

2250  T 1  =  E  K  P  ■:  -  X 
2260  T2=EXP<-U> 

2270  T 3  =  <  1  -T  1  > /  C  I-T2-) 

2  2 8 0  U< I >  =  D  m  < I >*X*T1/<P1*M2*A*< 1-T2>> 

2290  Ngmc < I ) =Ht ot * C T 1 -T2) /<. 1 -T2 > 

2 3 0 0  D 1  =  N tot* < X *T  1  *LOG < X > -T3*U*T2*L0G< U )  > ✓  P 1  ✓  C 1  - T 2  ') 

2310  Q  1  1  =  Q  1  1  +  D  1  2 
2  3  2  9  B  e  m  (  I  >  =  H  g  m  e  <  I  >  -  N  g  rn  c  f.  I  > 

2  3  3  0  S  q  d  =  S  q  d  +  D  e  m  <  I  >  '  2 
2340  Kl=Kl+Dl*Deo< I ) 

2350  NEXT  I 
2360  PRINT  “ SQD= " ; Sqd 
2370  C 1 =K 1 ,Q 1 1 
2  J  8  0  R  m  =  A  B  S  (  C 1  /  P 1 ) 

2390  IF  A  B  S  <  R  rn  >  <  =  .  2  THEN  2410 

2400  C  1  = .  2  R fu * C  1 

2410  IF  HBSCC1 )< . 001  THEN  2450 

2420  P 1 =  P 1 +  C 1 

2430  PRINT  "  P 1  = " ; P 1 

2440  GOTO  2200 

2 4 5 0  G 0 S U B  0 u t  put  2 

2460  DISP  “ P 1 =  ?  P 2  =  ?  BOTH  CONSTANT,  USING  MOTT  FORMULAS.  0,0  MEANS  SK IP"; 

2470  INPUT  Ml, M2 

2480  IF  M 1=0  THEN  GOTO  2660 

2490  Sqd=0 

25O0  PRINT  “  USE  CONSTANT  P1=";M1;"  AND  P2=";M2;"IN  MOTT  FORMULAS" 

251 0  V  i  1  >  =  M  u  p  <  1  )  2 
2520  D v ( 1 > = M u  p ( 1 > 

2  5  3 0  FOR  1=2  T  0  N  g 
2540  V (  I  >  =  M  C I > 

2550  D 11  (  I  '  =  D  m  (  I  > 

2560  NEXT  I 

2570  FOR  1  =  1  TO  N g 

2 5  8  0  X  =  <  V  <  I  >  /  M  2  >  <  1 ,  M 1  > 

2590  T 1 =  E  X P < - X ) 

2 6 0 0  W ( I >  =  D  v < I )*X*T 1  -  < M 1  * V < I > ) 

261 0  Ngrric  ( I  >  =  N t  ot  *T  1 
2620  IF  I < N g  THEN  2640 
2  6  3  0  N  g  rn  c  <  I )  =  N  ( I )  =  0 
2640  NEXT  I 
2650  G  0  S  U B  0  u  t  p  u t 2 

2660  DISP  "Pl=?  P2=?  ADJUST  P2  ONLY.  0,0  MEANS  SKIP"; 

2670  INPUT  M 1 , P 2 

2680  IF  M 1 =0  THEN  GOTO  2950 

2690  PRINT  “ADJUST  P2=";P2;"  ONLY,  WITH  CONSTANT  Pl=" ; M 1 

2700  Q  2  2  =  K  2  =  S  q  d  =  0 

2710  FOR  1=1  TO  Ng 

2720  X=  C  <  M  <  I  >  -Mm  i  n >  /  P2  >  <  1  /M 1  ) 

2  7  3  0  U  =  C  <  M  u  p  <  N  g  >  -  M  m  i  n  )  /  P  2  > A  <  1  /  M 1  > 

2748  T 1 =  E  X  P  < - X ) 

2750  T  2  =  EXP  < -U  > 

2760  T3=<  1-T1  1-T2) 

2770  W  <  I  >  =  D  m  (  I  >  *  X  *  T 1  /  <  M 1  *  (  M  (  I  >  -  M  m  i  n  >  *  (  1  -  T  2  )  ) 

2780  Ngmc  I  >  =Nt  ot  *  <  T 1 -T2  >  /  (  1-T2) 
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2790  D2  =  Ntot*CX*Tl-T3*U*T2)/Ml/P2/'C  1-T2) 

2.800  G.22»Q22+D2A2 

2810  DewC I >»NgmeC I >-Ngmc C I ) 

2820  Sqd=Sqd+DewC I >~2 
2330  K2=K2+DewC I)*D2 
2840  NEXT  I 
2850  PRINT  "SQD=" ; Sqd 
2860  C2=K2'Q22 

2870  Rm=ABSCC2/P2> 

2880  IF  Rm<=. 2  THEN  2900 

2890  C2= . 2/Rm*C2 

2900  IF  ABSCC2X.801  THEN  2940 

2910  P2=P2+C2 

2920  PRINT  "P2=“ ; P2 

2930  GOTO  2700 

2940  GOSUE  Output  2 

2950  STOP 

2960  END 

2970  Out  put  2:  Sqd=0 
2980  PRINT 

2990  PRINT  "  **NUMBER  WITH  MASS  GREATER  THAN  M  *  * " 

3000  PRINT  "  I  M< I >  EXPER  CALC  DEVIATION 

3010  FOR  1=1  TO  Ng 

3020  Dew=Ngmc  Cl)  -Ngrne  C  I  )  ■ 

3030  S  q  d  =  S  q  d  +  D  e u A  2 

3040  PR  I  NT  US  I NG  9@ ; I , M  Cl > , Ngme < I ) , Ngmc < I ) , Dew 

3O50  NEXT  I 

306O  PRINT  "SQD=";Sqd 

307O  DISP  "PRINT  NUMBER  IN  ITH  GROUP?  < 0 ,  1  >  "  ; 

3080  INPUT  I  rid 

3090  IF  I nd=0  THEN  GOTO  3210 

3100  PRINT 

3110  PRINT  "  **NUMBER  IN  ITH  GROUP** 

3120  PRINT  "  I  MCI)  EXPER  CALC  DEVIATION 

3130  S  q  d  =  0 

3140  FOR  1=1  TO  Ng 

3150  Nc  =W  C I >  *Nt  ot 

3160  Dew=Nc -N  C I > 

3170  Sqd=Sqd+Dev  ■  2 

3180  PRINT  USING  90; I , MC I > , NC I > , Nc , Dew 

3190  NEXT  I 

3200  PRINT  "800=" ; Sqd 

3210  PRINT 

3220  RETURN 
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